Collective dynamics of the high-energy proton-nucleus collisions 
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We analyze the proton-lead collisions at the LHC energy of ^snn = 5.02 TeV in the three-stage 
approach, previously used to successfully describe the relativistic A-A collisions. The approach 
consists of the early phase, modeled with the Glauber model, the event-by-event viscous 3-1-1 di- 
mensional (3-1-1 D) relativistic hydrodynamics, and the statistical hadronization at freeze-out. We 
show that features typical of collective dynamics, such as the harmonic flow and the ridge structures 
in the two-particle correlations in relative azimuth and pseudorapidity, may be naturally explained 
in our framework. In the proton-nucleus system the harmonic flow is generated from an initially 
event-by-event deformed system and is entirely due to these initial fluctuations. Notably, fluctu- 
ations of strength of the initial Glauber sources which yield the observed distribution of hadron 
multiplicities and, at the same time, lead to correct values of the elliptic flow coefficients both 
from the two- and four-particle cumulant method, as measured by the ATLAS collaboration. The 
azimuthally asymmetric flow is not modified significantly when changing the viscosity coefficient, 
the initial time for the collective expansion, or the initial size of the fireball. The results present 
an estimate of the collective component in the two-particle correlations measured experimentally. 
We demonstrate that the harmonic fiow coefficients can be experimentally measured with methods 
based on large rapidity gaps which reduce some of the other sources of correlations. 
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I. INTRODUCTION 

The recent interest in proton-nucleus (p-A) collisions 
stems from the expectations that the experimental data 
for this system could be used to test various theoret- 
ical approaches developed for relativistic collisions ll|, 
moreover, it could serve as a reference for experiments 
involving nucleus-nucleus (A-A) collisions. An interest- 
ing possibility is that the collective behavior clearly seen 
in the A-A collisions may be present already in the p-A 
collisions, and even in the proton-proton (p-p) collisions 
of highest multiplicity of the produced particles. The ex- 
perimental 0-Sl and theoretical p3l-[l3j investigations can 
provide a limit on the amount of the collective flovif in 
small systems, setting a boundary on the collective be- 
havior and helping to answer the questions: How small 
may the system be and under what conditions it is still 
describable with hydrodynamics? What observables are 
most sensitive to the collectivity? What is the inter- 
play of various stages of the dynamics, starting from the 
initial condition, through the intermediate evolution, to 
hadronization? 

The studies of the A-A collisions at RHIC jlJi] and the 
LHC |15| led to by now conventional interpretation of nu- 
merous observed phenomena in the mid-rapidity region 



via formation of a hot dense medium - the strongly inter- 
acting quark-gluon plasma - which evolves as a fluid and 
may be successfully described with relativistic hydrody- 
namics [16]. Basic phenomena supporting this view are 

• The harmonic flow (elliptic, triangular, higher- 
order) Bill, [13 

• Characteristic ridge structures seen [J, [l^ [I^ in 
the 2-particle correlation functions in relative pseu- 
dorapidity and azimuth, possible to explain with 
harmonic flow [ll,[23|. 



Specific features of the interferometric radii |21|.|2 



Jet quenching 
therein) . 



(see, e.g., Q and references 
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In the study presented in this article we investigate the 
first two items from the above list for the case of the 
p-Pb collisions, recently studied at the LHC at the colli- 
sion energy of y^J/vw = 5.02 TeV. We apply a treatment 
based on hydrodynamics to find quantitative estimates 
for the measured quantities, extending an early analysis 
of flow Tl] by one of us (PB) in the p-A and deuteron- 
nucleus systems, as well as the more recent event-by- 
event studies of the correlation functions [IJ] and the 
interferometric radii [24| . As will turn out, our results 
agree at a quantitative or semi-quantitative level with 
the experimental data for the highest centrality classes, 
supporting the collective picture of the most central p-Pb 



collisions. Our results also set the background for more 
elementary explanations of the correlation studies, based 
on saturation and the color-glass-condensate (CGC) the- 
ory |25l430| . We note that certain features can also be 
obtained with cascade models [3l| . 

We note that a certain degree of collectivity has been 
suggested for the p-p collisions as well [y, il(>|, where a 
same-side ridge is observed in the 2D correlations func- 
tions 8] for the highest multiplicity events. This may in- 
dicate the presence of azimuthal correlations in the gluon 
emission from the initial state [25, 26, 28 30]. However, 
the same-side ridge observed in the p-p collisions could 
also result from the collective expansion of the created 
medium [10| . The intriguing questions concerning the p- 
p collisions will not be explored in this work, devoted to 
the detailed analysis of the p-Pb case. 

Finally, we stress that our method is applicable to soft 
physics, related to particles produced with transverse mo- 
menta lower than, say, 2 GeV. 



II. THE THREE-STAGE APPROACH 

Our event-by-event approach is based on, by now, a 
standard picture involving three stages: generation of the 
initial densities, hydrodynamic evolution, and hadroniza- 
tion. Certainly, variants of the modeling of each stage are 
present in the literature. We use the Glauber approach 
as implemented in GLISSANDO [33] to model the ini- 
tial phase, the 3-fl D viscous hydrodynamics [33], and 
the statistical hadronization as implemented in THER- 
MINATOR M. 



A. Initial conditions 

The initial condition is generated with GLIS- 
SANDO i32], implementing various variants of the 
Glauber model [33, [Sa]. The parameters of the calcu- 
lations are similar as in [ll|, except that they are ad- 
justed for the collisions energy of i/sjvjv = 5.02 TeV. 
Thus we take the nucleon-nucleon (NN) inelastic cross 
section a^N = 67.7 mb, moreover, we use a realistic 
(Gaussian) wounding profile [33] for the NN collisions. 

In the Glauber model, when a NN collision occurs, a 
source is produced, meaning a deposition of energy in 
a location in the transverse plane and spatial rapidity. 
In the standard wounded nucleon model it is assumed 
that a point-like source is located at the center of each 
participating nucleon, which leads to a rather large initial 
transverse size in the p-Pb system. Locating the source 
in the center-of-mass of the NN pair is also a possible 
model choice; it leads to smaller initial distributions. In 
the results presented below we use both variants of the 
model, termed standard and compact. That way we may 
estimate the uncertainty related to modeling the initial 
phase within the Glauber treatment. 



TABLE I. The mean values and standard deviations of the ba- 
sic characteristics of the initial distributions for the centrality 
class 0-3.4%. 





standard 


compact 


Glauber-fNB 




mean std. dev. 


mean std. dev. 


mean std. dev. 


{rY^^ [fm] 


1.54 0.15 


0.93 0.06 


1.45 0.22 


£2 


0.25 0.13 


0.19 0.09 


0.34 0.16 


es 


0.29 0.13 


0.18 0.09 


0.32 0.15 



There is another important effect in the initial stage 
that influences the results: The weight of each source 
fluctuates according to some statistical distribution, sim- 
ply reflecting the fact that each NN collision may produce 
a different number of partons and therefore lead to vary- 
ing deposition of the entropy. In the simulations 32] , this 
feature is achieved by overlaying a suitable distribution 
of strength w over the spatial distribution of the partic- 
ipant nucleons. As described in Sect. Ill El the observed 
multiplicity distributions can be described as convolution 
of the number of participant nucleons and a negative- 
binomial distribution. At the stage of the formation of 
the initial fireball it is equivalent to imposing fluctua- 
tions of the entropy deposited per participant nucleon 
following the F distribution. 



PrH 



Tin) 



(1) 



with K — 0.9, as it leads to correct multiplicity distribu- 
tion of the produced particles, cf. Sec. Ill El This case is 
labeled Glauber+NB as it eventually gives the multiplic- 
ity distribution as a convolution of the Glauber Monte- 
Carlo distribution of participant nucleons and a negative 
binomial distribution. It is the most physical one and 
leads to best results when compared to the experiment. 
By construction, in the Glauber -l-NB case we place the 
sources at the centers of the participant nucleons, as in 
the standard case. 

In all cases, after generating the spatial positions of 
sources, we smooth them with a Gaussian profile of width 
0.4 fm. This physical effect (the sources do have a non- 
zero width) is also essential for hydrodynamics, which 
requires sufficient smoothness of the initial conditions. 
The smoothed initial distribution is placed on the 3+1 D 
lattice with spacing of 0.15 fm and then the event- by- 
event hydrodynamics is run. 

Some features of the resulting initial distributions are 
shown in Figs. [T] and [5] for the collisions at the high cen- 
trality, c=0-3.4%. We use here a few hundred configu- 
rations generated with GLISSANDO which are later fed 
into the event-by-event hydrodynamic evolution. The ba- 
sic properties of the distributions are listed in Table. HI 




FIG. 1. (Color online) The distribution of the transverse rms 
radius for the initial configuration in the standard case (solid 
line), compact case (dashed line), and the Glauber+NB case 
(dotted line), for the centrality class 0-3.4%. 
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FIG. 2. (Color online) The distributions of the scaled ec- 
centricities e„ and scaled harmonic flow coefficients «„{2}, 
n = 2, 3 for the standard case. The thin solid line shows the 
Wigner distribution of Eq. Q. The flow coefficients Vn{2} 
are discussed in the following sections. 



B. Initial size 

At first glance, a rather surprising feature is the large 
transverse rms size of the initial distributions. To under- 
stand, consider first the standard case, where the sources 
are located in the centers of the proton and of each of 
the participants from the lead nucleus. If the geometric 



"hard-disk" of radius R were used for wounding, then 
the inelastic cross section would be a^N = ttR^, from 
where R = 1.47 fm, corresponding to rms of 0.98 fm in a 
single NN collision. This would be the uncertainty of the 
location of the point-like source created in the NN colli- 
sion in the model. However, we use the realistic Gaussian 
wounding profile [37| of the form 



p{b)^Aexp{-7TAb^/aNN), ^ = 0.92. 



(2) 



After folding with the distributions of the nucleons in 
the Pb nucleus (via the Monte Carlo procedure in GLIS- 
SANDO) and after smoothing the positions with Gaus- 
sians of width 0.4 fm located in the centers of each source, 
we obtain the rms radius of the initial distribution listed 
in Table [H namely 1.54 fm. 

In the compact case, where the the sources are placed 
in the centcr-of-mass of the colliding proton and the nu- 
cleon from the Pb nucleus, the source is significantly 
smaller, with the transverse rms of 0.93 fm. More in- 
volved models of the initial stage have been consid- 
ered |3g| . The details of the energy deposition, such as 
the fluctuations and small scale structures, are relatively 
more important in p-Pb than in A-A collisions. Our cal- 
culation, which uses two cases with significantly different 
initial size of the fireball can serve as an illustration of the 
effects of the variation in initial size on the final harmonic 
flow observables. 



C. Initial eccentricities 

The participant eccentricities shown in Figs. [5] are de- 
fined in the usual way for a given event as 



J dxdyr"'p{x, y)e^'^ 
J dxdyr^p{x,y) 



(3) 



where p{x, y) is the initial transverse entropy distribu- 
tion in the fireball at zero pseudorapidity and V^n is the 
azimuthal angle of the event plane. For the p-Pb system, 
the origin of the non-zero eccentricity lies in the fluctua- 
tions of the positions of the participant nucleons. ^From 
the formulas of Appendix D of Ref. 36] it is straight- 
forward to obtain the result that for the very central 
collisions, c = 0, where the average distribution in the 
Pb nucleus seen by the proton is azimuthally symmetric, 
the scaled eccentricities s — e„/(e„) calculated from the 
positions of the participant nucleons follow the Wigner 
distribution 



w{s) 



TTS 



■exp( 



TTS 



(4) 



independently of the rank n of the harmonic component. 
The distribution has (s) — 1 and var(s) = 4/7r — 1 [36|. 
We note that this universality is clearly seen in Figs. [U 
up to the statistical noise and a slight departure from 
the central case of c = 0. The eccentricities are calcu- 
lated from the smooth coarse-grained lattice distributions 



which introduces some small corrections with respect to 
the eccentricities calculated from discrete positions of the 
the participant nucleons |39l |. 



D. Pseudorapidity distribution 

It is assumed that the initial transverse and longitudi- 
nal distributions are factorizcd. This assumption plays 
a key role in the interpretation of the development of 
the ridge structures in the hydrodynamic approach. It 
means that the transverse distribution is, within a rea- 
sonable range around the central region, independent of 
the pseudorapidity, i.e., approximately boost invariant. 
This leads to a correlation of "geometry" for the fire- 
ball slices separated by A77, and, in consequence, to the 
correlation of flow. If an (approximately) boost-invariant 
fireball is formed, azimuthal correlations due to collective 
flow show up [l7l.[T9|. 

For the shape of the longitudinal distributions in the 
NN center of mass frame we use the following profiles in 
the space-time rapidity Tyn : 



/(r/||)± =exp 

iVb ± V\ 



(l^|||-%)^ 



2^2 



h\\-m)j X 
'-d{yb±v\\), (5) 



Vb 

with rjQ — 2.5, cr,, = 1.4, and yh = 8.58 denoting the 
beam rapidity. The indices + and - correspond, respec- 
tively, to the distribution generated by the forward and 
backward moving participant nucleons. The same func- 
tional form of the profile has been successfully used in 
Refs. [40J to describe features of the A- A collisions, in 
particular the spectra in pseudorapidity and the behav- 
ior of the directed flow at RHIC. A phenomenological 
motivation for such "triangular" parameterizations has 
been discussed in [40, |41| . The resulting long-range cor- 
relations in pseudorapidity are strong in the asymmet- 
ric p-Pb collisions, hence the reaction planes at different 
pseudorapidities are aligned (Sec. IIVB|) . To summarize, 
the initial conditions for hydrodynamics are the product 
of the smoothed transverse Glauber distribution in the 
transverse plane and the function ([5|) in the longitudinal 
direction. 



E. Multiplicity distribution and fireball 
fluctuations 



The Glauber Monte Carlo approach provides event- 
by-event fluctuations in the number of NN collisions and 
their distribution in the transverse plane. This mecha- 
nism explains most of the observed fluctuations in the 
shape of the fireball in the A- A reactions [i^, |43| . The 
event-by-event hydrodynamic expansion of the fluctuat- 
ing fireball generates azimuthally asymmetric flow and 
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FIG. 3. (Color online) Multiplicity distribution of tracks witii 
Pi_ > 0.4 GeV and |J7| < 2.4 measured by CMS g^]. The dot- 
ted and solid lines denote the convolution of the distribution 
of participant nucleons from GLISSANDO with the Poisson 
and negative binomial distributions, respectively. 

its fluctuations [13, [H, IH, \I^ . Some observables indi- 
cate that additional sources of fluctuations are present, 
beyond the fluctuations in the number of participant 
nucleons, e.g., the event-by-event distribution of har- 
monic flow coefficients [15, 45j| or the multiplicity dis- 
tributions |38l.l46j. 

It is well known that the multiplicity distributions in 
the p-p collisions can be described by the negative bino- 
mial distribution 14811 



NxAn) 



r(n + K)A"K'' 
F(K)n!(A + K)"+''^ 



(6) 



where the multiplicity n has the mean and variance given 
by A and A(l -I- X/k), respectively. Below we argue that 
in the p-Pb collisions we have a similar situation. In 
Sec. Ill Al we stated that one should overlay a weight dis- 
tribution over the spatial distribution of participant nu- 
cleons - a feature implemented in GLISSANDO. Now we 
show how this distribution can be adjusted in such a way 
that the multiplicity data are properly reproduced. 

We use the multiplicity distributions of tracks observed 
in the minimum-bias p-Pb collisions !47i]. In the three 
stage model of particle production, multiplicity fluctua- 
tions come from the fluctuations of the initial entropy of 
the fireball, from the entropy production during the vis- 
cous hydrodynamic stage, and from the statistical emis- 
sion of hadrons at freeze-out. To a good approximation, 
in the considered regime of centralities the entropy after 
the hydrodynamic expansion is directly proportional to 
the initial entropy, which refiects the fact that the de- 
terministic hydrodynamic evolution does not introduce 
fluctuations [4^. For independent statistical emission, 
the number of emitted hadrons is proportional to the fi- 
nal entropy and follows the Poisson distribution pM| . 

First, we consider the case where there is no weight 
distribution overlaid over the participant nucleons. Then 
the initial entropy is proportional to the number of par- 
ticipant nucleons and the distribution of the observed 
tracks is a convolution of three distributions: the dis- 



tribution of participant nucleons, a Poisson distribution 
for each participant with a mean A defined as the aver- 
age number of particles produced per participant, and 
a binomial distribution with success rate p giving the 
probability of recording a track in the CMS acceptance. 
The folding yields the multiplicity distribution of the pro- 
duced hadrons of the form 



pPb 5020GeV A/part=19 



P(n)=^Pp,rt(z) 



(Api)" e~^P' 



(7) 



where Ppartii) is the distribution of the participant nu- 
cleons from the Glauber Monte Carlo. The parameter 
Ap — 5.36 is chosen to reproduce the mean number of 
the observed tracks. As we can see (the dotted line in 
Fig. [3]), the multiplicity distribution from the Glauber 
model convoluted with the Poisson distribution is much 
too narrow and does not reproduce the experimentally 
observed high-multiplicity tail. 

The above shows that inserting a distribution of 
weights over the participant nucleons is necessary. In 
that case the distribution of the observed tracks is a con- 
volution of four distributions: the distribution of par- 
ticipant nucleons, an overlaid distribution of weights, a 
Poisson distribution for the production of hadrons, and 
a binomial distribution for the experimental acceptance. 
When we use the distribution of weight in the form of 
the r distribution ([T]), then its folding with the Poisson 
distribution yields the negative-binomial distribution ([6]) , 
which we now take for hadrons produced per participant 
nucleon. One finds 



P{n) = y^,Ppart{i)N'pXt.Ki{n). 



(8) 



The experimental multiplicity distribution is now very 
well reproduced with the parameter values Xp — 5.36 and 
K = 0.9 (sohd line in Fig. [3]). We refer to this calculation 
as Glauber+NB. 

The procedure outlined above is clearly a simplified 
picture of the multiplicity fiuctuations in the relativis- 
tic nuclear collision. Further effects could be important, 
in particular the shape of the multiplicity distribution 
can depend on the pseudorapidity window, the track ac- 
ceptance of the CMS detector is not uniform, particles 
from jets that contribute to the the tails of the multiplic- 
ity distribution do not increase the fiuctuations in the 
thermalized fireball, or the entropy production in vis- 
cous hydrodynamics is not exactly linear. Nevertheless, 
the considered mechanism of additional density fiuctua- 
tions in the fireball can serve as a model to illustrate its 
expected effects on the eccentricities and the harmonic 
flow coefficients. Such effects in A-A collisions have been 
considered previously in the Glauber scheme [36l |50| and 
found to be relevant. 
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FIG. 4. (Color online) A sample event evolution, visualized 
via the freeze-out isotherms in the x — t plane (solid) and the 
y — T plane (dashed). The standard case. 



pPb5020GeV A/part=19 



o 
E 4 




1 



2 4 

X (y) [fm] 

FIG. 5. (Color online) Same as Fig.[4]for the compact-source 
case. 



F. Hydrodynamics 

The hydrodynamic model used in this work is de- 
scribed in detail in [llj. It carries out an event- by-event 
3-|-l D evolution and includes the shear and bulk viscosi- 
ties. The multiplicity expected in central p-Pb collisions 
is extrapolated linearly in the number of the participating 
nucleons from the minimum bias results of the ALICE 
collaboration [4]. That way the average initial entropy 
per participant is adjusted. The shape of the entropy 
distribution follows the distribution of sources described 
in the previous section. The starting time of hydrody- 
namics is chosen to be tq = 0.6 fm/c for the standard 
case, but the evolution with the choice tq = 0.2 fm/c 
is also studied. The relatively short total duration of 
the collective expansion phase makes the results more 
sensitive to this very early stage a nd, possibly, to some 
nonequilibrium transient behavior [51| . The ratio of the 
shear viscosity r] to entropy density s is 77/s = 0.08 or 
77/5 = 0.16, while the ratio of the bulk viscosity C to 
s in the hadronic phase is C,/s = 0.04 [52|. For each 
case (standard, compact, higher viscosity, lower initial 
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FIG. 6. (Color online) The model predictions for the mid- 
rapidity transverse momentum spectra of tt"^ for the most 
central events. The solid line is for the standard calculation, 
with the initial source rms size of 1.54 fm, r]/ s = 0.08, and 
the initial time tq = 0.6 fm/c. The dotted line shows the re- 
sults for the initial time of 0.2 fm/c (and the other parameters 
unchanged), the dashed line stands for the calculations with 
•q/s = 0.16, and the dash-dotted line represents the calcula- 
tion with the initial rms size 0.93 fm/c. The solid line with 
the triangle symbols shows the Glauber-|-NB results, the case 
where the Glauber Monte Carlo initial conditions are convo- 
luted with the r distribution (cf. Sec. HIE) . 



time, or Glauber -|-NB) we produce initial configurations 
that are evolved event-by-event with hydrodynamics to 
obtain freeze-out hypersurfaces of the constant tempera- 
ture T/ = 150 MeV. 

Two typical evolution histories for the standard and 
compact case, both for A^part = 19, are depicted in 
Figs. I311S1 where we show isotherms at Tf = 150 MeV in 
the x — T plane (solid lines) and in the y — r plane (dashed 
lines). We note that although the systems originally have 
difi^erent sizes, the spatial spread of the isotherms at later 
times is similar, about 5 fm. The evolution of the stan- 
dard source is about 15% longer than from the compact- 
source case. The radial flow is larger for the compact 
case, as the system is more squeezed initially. This leads 
to 20% higher values of the average transverse momen- 
tum. 

The transverse momentum spectra depend on the 
choice of initial conditions. The pr-spectra for tt^ are 
shown in Fig. [B) We notice that the spectra get harder 
with the increase of the shear viscosity, the decrease of 
the initial time tq , or in the Glauber -l-NB case. The hard- 
ening of the spectra is most pronounced when starting 
the calculation from the compact source, which involves 
larger gradients present in the system. 



G. Statistical hadronization 

For each freeze-out configuration we generate 1000 
THERMINATOR events to efliciently improve the statis- 
tics. This is a technical point, as physically one event 
should hadronize into one set of hadron distribution. The 
trick of running multiple THERMINATOR simulations 
on the same hydro event allows us to efliciently improve 
the statistics, as the computing time for the hydro evo- 
lution is two orders of magnitude longer than for the 
generation of a single THERMINATOR event. 

We note that the statistical hadronization built in 
THERMINATOR contains the non-fiow correlations 
from all resonance decays. The use of a full Monte Carlo 
generator of hadron distributions is also of practical mer- 
its, as it allows implementation of the kinematic cuts, ac- 
ceptance or efficiency from the experimental setup, which 
is crucial in comparisons to the data. 



H. Local charge conservation 

Sizable correlations among opposite-charge particles 
result from the local charge conservation [5J]. There 
are indications that this eflect is generated at hadroniza- 
tion ^, _54J, i.e., at the late stage of the reaction. Our 
implementation of the charge balancing is based on the 
assumption that the particle- antiparticle pairs of charged 
hadrons are emitted locally at freeze-out, carrying ther- 
mal distribution. The mechanism is described in detail 
in Ref. [2fl]. 



I. Transverse-momentum conservation 

In our studies of the correlation variables we enforce 
the global transverse momentum conservation, which is 
important in correlation analyses ^55,] . In particular, it af- 
fects the shape of the two-particle correlations in relative 
pseudorapidity and azimuth. To satisfy the constraint 
approximately we require the following condition on the 
global transverse momentum: 



\ 



E^ 



E 



Py 



<Pt 



(9) 



where i labels particles in the event. We have found nu- 
merically that in the central p-Pb system it suffices to 
take Pt = 5 - 10 GeV. That way we retain 5-10% of the 
least-Pr events from our full sample. A further lowering 
of Px does not affect the correlation results, while it de- 
teriorates the statistics. The momentum conservation is 
imposed when calculating the di-hadron correlation func- 
tions. For the calculation of the elliptic and triangular 
flow, imposing the momentum conservation in that form 
is irrelevant. 



J. Centrality definition 



The simplest determination of the centrahty classes 
in our model can be obtained from conditions on the 
number of the participant nucleons. The collisions with 



N, 



part 



> 18 amount to 3.4% of the most central events 



from GLISSANDO. The next most central class is defined 
as 16 < Npart < 17, which forms centrality 3.4 — 7.8%. 
On the other hand, as noted in [ll|, simplistic centrality 
definitions based on the impact parameter are ill defined 
for central p-Pb collisions. 

For the Glauber+NB case we define the centrality 
classes not by the number of participant nucleons, but 
through the total initial entropy in the fireball, i.e., we 
take into account the fluctuations from the overlaid F 
distribution. A more accurate determination, following 
closely to the experimental setup, should impose the cuts 
on the flnal multiplicity of the produced particles, instead 
of on Npart or the initial entropy. 



III. TWO-PARTICLE CORRELATIONS 
A. Definitions 

The basic objects of the study of this section are the 
two-dimensional two-particle correlation functions in rel- 
ative pseudorapidity and azimuth. These quantities in 
comparison to the CMS data [J| have already been ana- 
lyzed in 12]. Here we extend this analysis, comparing to 
the ATLAS data H as well. 

The simplest definition of the correlation function in 
the considered kinematic variables is 



C(A,7,A(/.) 






d^Afpo 



(10) 



\dAjjdA0/ mixed events 



If the correlations were absent, C(Ar/, A(/)) = 1, thus 
unity is a natural scale for this measure. This correlation 
is used by the ATLAS collaboration ^Q]. 

The "per trigger" correlation function, used by the 
CMS collaboration, is defined as B 



Ctrig(A?7,A<?i) 



1 d^NP""" 
N dAr] dA(t) 



5(0,0) 



g(A77,A^) 
B(A?7,A</>)' 

(11) 



with A?7 and Acf) denoting the relative pseudorapidity and 
azimuth of the particles in the pair. The signal and the 
mixed-event background are defined with the pairs from 
the same event, and the pairs from the mixed events, 
respectively: 



S{Af],A4)) : 

B(Ar/,A0) 



, 1 d'^NP''" 
'NdATjdAc/)' 
1 d'^NP^" , 
^N dAridAcI)- 



(12) 



c=0-3.4%, 0.5<p^<4 GeV 




FIG. 7. (Color online) The charge-independent correlation 
function C(Ar;, A(j!>). The total transverse momentum is ap- 
proximately conserved with the condition Pt < 10 GeV. 
Charge balancing is included. 



The number of particles N in the prefactor (denoted by 
CMS as A^trig) is the number of charged particles in a 
given centrality class and acceptance bin, corrected for 
the detector efficiency. The multiplication with the cen- 
tral bin content B{0,Q) in Eq. (jlip brings in the inter- 
pretation of Eq. (jlip as the average number of correlated 
pairs per trigger particle. 

To make quantitative comparisons easier, one also uses 
the projected correlation functions. A function used by 
the ATLAS collaboration is defined as 



Y{A^) 



J B{A(j))d{A(t)) 



ttN 



C{A^)-bzYAM, (13) 



'mixed events- 



where S'(A0) and i3(A0) are averages of 5(Ary, A0) and 
B{Ari, Ac/)) over the chosen range in Ary avoiding the cen- 
tral region, in particular 2 < \Ari\ < 5 in the ATLAS 
analysis, and the constant ^zyam is such that the mini- 
mum of Y(A(f>) is at zero. 



B. Comparison to the ATLAS data 

The result of our simulations for the most central p-Pb 
collisions (c = — 3.4 %) with the kinematic cuts corre- 
sponding to the ATLAS setup Q is shown in Fig. [T] We 
display the standard-source case, as for the compact or 
Glauber -l-NB cases the results are quantitatively similar. 
We note the two prominent ridges, generated with flow, 
as well as the central peak, coming in our simulation from 
the charge balancing |2(j |. 

The same-side ridge appears naturally as a conse- 
quence of the collective flow. More precisely, in our 
framework the shape and flow in the fireball in the for- 
ward and backward rapidity regions is correlated, reflect- 
ing the assumption on the factorization of the transverse 
and longitudinal distributions in the initial condition. In 
particular, the principal axes of the elliptic flow are cor- 
related along the whole pseudorapidity span. Thus, there 




c=0-3.4% 

0.5 < p^ < 4.0 GeV 
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FIG. 8. (Color online) Projected and ZYAM-subtracted 
correlation function YA{(f)) for the most central p-Pb colli- 
sions for the standard source (solid line) and compact source 
(dashed line), compared to the ATLAS data (points) at 
T,E^'' > 80 GeV. The total transverse momentum is approx- 
imately conserved with the condition Pt < 5 GeV. Charge 
balancing is imposed. 
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FIG. 9. (Color online) Projected correlation function CA{(j}) 
for the most central p-Pb collisions for the standard source 
and several values of the maximum total transverse momen- 
tum Pt, listed in the legend. 



Effects of transverse momentum conservation 



are more pairs with A0 ^ and A0 ^ tt regardless of 
A77. This "flow explanation" of the ridge formation is 
appealing in its simplicity. 

Next, to compare quantitatively to the data, we look 
at the projected correlation function Y(A(f)). There is a 
technical issue which must be discussed. By construction, 
the prefactor oiY{A(j)) is proportional to {N{N—1))/{N) 
- the ratio of the average number of pairs to the aver- 
age number of particles. Thus to reproduce Y{A(j)) in a 
model calculation one needs to have proper correlations, 
but also correct fluctuations of the multiplicity. The sec- 
ond requirement is not easy to accomplish, in particular 
as ATLAS is using the transverse energy to define the 
centrality classes, and the fiuctuations of multiplicity are 
large. For that reason in our comparison we rescale our 
model Y{A(j)) in such a way that the subtraction con- 
stant 6zYAM is the same in the model and in the experi- 
ment. This assumes that the mechanism generating the 
fiow and the ridge structures is "factorisable" from the 
multiplicity fluctuations. 

The result of this procedure is shown in Fig. [5] for the 
most central collisions. We note that experimental data 
fall within our model results for the standard (solid line) 
and compact (dashed line) sources. We note that the 
compact source, leading to larger flow, has more promi- 
nent ridges. 

The CMS correlation data for the p-Pb collisions have 
been compared to in our previous short paper |24| , hence 
we do not repeat these results here, but only mention 
they are in semiquantitative agreement with the data. 



We can now demonstrate the relevance of the trans- 
verse momentum conservation and the simple procedure 
introduced in Sec. Ill II We use the projected correlation 
function C(A0) for that purpose. We note that limiting 
the value of the maximum total transverse momentum 
Pt in the accepted events moves the strength from the 
same-side ridge to the away-side ridge. This is natural, as 
the momentum conservation increases the back-to-back 
motion of the particle. We note that for a practical pur- 
pose it is enough to use Pt < 5 — 10 GeV. A further 
reduction changes the results very little at the expense 
of deteriorating the statistics. The numerical results, dis- 
playing the mentioned convergence, are shown in Fig. |9l 



IV. HARMONIC FLOW 
A. Cumulant method 

^From the two- and four-particle cumulant method [56| 
we obtain the values of the flow coefficients collected 
in Table [TTl The kinematic cuts correspond to the AT- 
LAS experimental setup. We compare the standard and 
the Glauber -l-NB simulation, without or with the local 
charge conservation, and list the results for the two high- 
est centrality classes. In Table IIIII we show the depen- 
dence of W2{2} for the most central events on the param- 
eters of the model: the value of the sheer viscosity rj and 
the time when hydrodynamics is initiated, tq. 

We note several qualitative features from Tables [III and 

mi 

• The dependence on centrality is very weak, as ex- 
pected from the flow generated mainly by the fiuc- 
tuations of the initial condition. 



TABLE II. Model predictions for the elliptic and triangu- 
lar flow coefficients from the cumulant method for the p- 
Pb collisions at ^sJTiv = 5.02 TeV. The cuts \ri\ < 2.5, 
0.3 < pt < 5 GeV correspond to the ATLAS setup. The stan- 
dard and Glauber-j-NB cases are used, with 77/s = 0.08 and 
To = 0.6 fm/c, without and with charge balancing. The errors 
are statistical and reflect the accumulated number of simu- 
lated THERMINATOR events (the actual error are somewhat 
larger due to a small number of the event sample). 



std., c=0-3.4%, 0.3<p^<5GeV 





c=0-3.4% c=3.4-7.8% 


standard, no balancing 


V2{2y [10-^] 


3.70(1) 3.78(2) 


«3{2}^ [10-3] 


1.04(1) 0.95(1) 


«2{4}* [lO-"^] 


-0.4(4) 1.83(5) 


«3{4}^ [10-«] 


0.0(2) -0.3(3) 


Glaubcr+NB, no balancing 


V2{2r [10-3] 


8.18(12) 8.24(10) 


V3{2r [10-3] 


1.52(8) 1.51(6) 


V2{4y [10-«] 


15(7) 16(6) 


«3{4}* [10-«] 


-2(2) -2(2) 


Glauber-j-NB, with balancing 


V2{2r [10-3] 


8.22(7) 8.68(6) 


«3{2}^ [10-3] 


1.57(4) 1.62(4) 


«2{4}* [10-«] 


19(4) 19(4) 


«3{4}* [10-S] 


-1(1) 0(1) 



TABLE III. Parameter dependence of the predictions for 
the elliptic and triangular flow coefficient from the two- 
particle cumulant method for the p-Pb collisions at y/SNN ~ 
5.02 TeV, c = - 3.4 %, |?7| < 2.5, 0.3 < pt < 5 GeV. Charge 
balancing not included. 





V2{2} [%] «3{2} [%] 


standard 


77/s = 0.08, TO = 0.6 fm/c 
77/s = 0.08, TO = 0.2 fm/c 
77/s = 0.16, TO = 0.6 fm/c 
77/s = 0.16, TO = 0.2 fm/c 


6.09(1) 3.22(2) 
7.44(1) 4.49(1) 
5.57(1) 2.67(2) 
7.12(2) 4.01(2) 


Glauber-hNB 


77/s = 0.08, TO = 0.6 fm/c 


9.0(1) 3.9(2) 



• The ■(;„{4}'' coefficients are, within the statistical 
limit of our simulations, compatible with zero for 
the standard case, while for the Glauber-j-NB sim- 
ulations t'2{4}'* is positive. This again shows the 
fluctuation nature of the generated flow from the 
Glauber initial conditions. Additional fluctuations 
of the entropy deposited initially in the fireball in- 
crease the eccentricity and yield a nonzero value of 
W2{4} (cf. Table U). 

• Increased sheer viscosity quenches, as expected, the 
flow. The relative effect is stronger for higher har- 
monics. 




FIG. 10. (Color online) The flow coefficients 112 (A77) (up- 
per lines) and V3{Arj) (lower lines) calculated from the two- 
particle correlations as function of the relative pseudorapidity 
of the particles in the pair. The solid and dashed lines are for 
the unlike- and like-sign pairs, respectively. The central peak 
is due to charge balancing and, to a lesser extent, resonance 
decays. 



• A shorter time of starting hydrodynamics increases 
the flow, which again is expected. 

• The effect of the local charge balancing increases 
somewhat the flow coefficients. 

As a matter of fact, the first two items above are cru- 
cial for the proper interpretation of the observed phe- 
nomenon. Detailed comparisons of the model predictions 
to experimental measurements provide a way of learning 
about the shape and fluctuations of the initial density in 
the p-Pb system. The observation of nonzero ■y2{4} by 
the ATLAS Collaboration 57[ indicates that in the small 
interaction region formed in the p-Pb collisions the large 
fluctuations of the energy deposited in each NN collision, 
as present in the Glauber-j-NB case, are crucial. Thus 
the initial conditions from the Glauber model in p-Pb 
collisions are fluctuation-dominated, analogously to the 
central A-A collisions. The same observation applies to 
the final elliptic and triangular flow in p-Pb collisions. 

In view of the recent experimental results for the 2D 
correlations functions in p-Pb collisions, it is interesting 
to look at the possibility of measuring directly the har- 
monic flow coefficients. We plot the elliptic and triangu- 
lar fiow coefficients as functions of the pseudorapidity gap 
in Fig. [TOl The quantities are obtained in our hydrody- 
namic model from the Fourier decomposition of the cor- 
relation function Ctrig(A77, A0) [20]. The non- fiow effects 
present in our model are important only for pairs of small 
pseudorapidity separation. In the intervals |A?7| > 2 the 
non-flow effects from the resonance decays and the local 
charge conservation can be neglected. We note that the 
flow coefficients in Fig. [TU] are sizable, thus could be mea- 
sured. It must be noted, however, that other sources of 
non-flow correlations may be present also in that kine- 
matic region, but with smaller amplitudes, as measured 
in the p-p collisions Q. 
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FIG. 11. (Color online) The p± dependence of the elliptic flow 
coefficient of charged particles with \ri\ < 1, obtained from the 
second-order cumulant method. The solid line corresponds to 
the standard calculation, rj/s = 0.08, ro — 0.6 fm/c, with 
the initial source rms size of 1.5 fm, the dotted line shows 
the case where the initial time is reduced to 0.2 fm/c, the 
dashed line stands for the calculations with increased viscos- 
ity, r)/s = 0.16, and the dash-dotted line represents the case 
of the compact source with rms size of 0.9 fm/c. Finally, the 
solid line with the triangle symbols shows the Glauber-|-NB 
case (Sect. IIIE|I . 



The pj^-dependent elliptic and triangular flow coefR- 
cients calculated with the two-particle cumulant method 
|56| are presented in Figs. [TT] and [T^ In the px < 2 GeV 
range, where hydrodynamics applies, the flow coefficients 
show a typical hydrodynamic behavior and the magni- 
tude of the flow is large. We find the elliptic (triangular) 
flow of about 10% (5%) for pj_ - 1 GeV. The resuhs 
are sensitive to the physical parameters of the model 
(cf. Table IIII[) . The flow decreases for larger viscosity 
or when using compact initial conditions, and increases 
when starting the hydrodynamic evolution earlier. It also 
increases with the presence of additional initial fluctua- 
tions, as in the Glauber -l-NB case. We notice a larger 
relative variation for the triangular flow than for the el- 
liptic flow when varying the parameters. 



FIG. 12. (Color online) Same as Fig. [TT] for the triangular 
flow coefficient, W3{2}. 
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FIG. 13. (Color online) The elliptic flow coefficient of charged 
particles for \ri\ < 2.5, 0.3 < p± < 5.0 GeV from the cumu- 
lant method ?;2{2} and i'2{4}, and from the di-hadron corre- 
lation function measured by the ATLAS Collaboration [53], 
compared to our hydrodynamic calculation for the standard 
case (u2{2} at centralities 0-3.4% and 3.4-7.8%) and for the 
Glauber-fNB case (^2(2} and W2{4} at centralities 0-5%, 5- 
10% and 10-20%). The corresponding transverse-energy for 
the centralities in the model calculations is obtained via in- 
terpolation of the experimental values. 



B. Scalar product method 

The correlation between particles produced in p-Pb 
collisions can have different origin. A way to reduce some 
of the non-flow contributions to the harmonic flow coeffi- 
cients is to use methods involving a rapidity gap between 
the reference particles defining the event plane orienta- 
tion and the particles used to calculate the flow coeffi- 
cient. In this subsection we present results for the scalar 
product method [53, [5^1 • One defines the Q„ vector 



Q 



A,B i-i/n 



E 



WkC 



incfik 



(14) 



as a sum over charged particles in a given reference bin 
{A or B). We use two definitions of the event plane, 
one with charged particles with 0.3 < p± < 3 GeV and 
2.0 < ?7 < 2.5 (Pb side), or with -2.5 < 77 < 2.0 (proton 
side). The weights are equal to the transverse energy 
{wk = E±) for the 3.2 < \ri\ < 4.8 bin. The resolution 
correction is 



Qi. 



'{Q^QS){Q^Q^) 



{Q 






(15) 



where the reference bin C is defined in all cases as 0.3 < 
p± < 3 GeV and \ri\ < 0.5. We have checked that the 
results do not differ noticeably when changing the p±^ or 
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FIG. 14. (Color online) The elliptic (solid symbols and dashed 
line) and triangular (open symbols and solid line) flow coeffi- 
cients obtained with the scalar product method. The circles 
and squares represent the calculation using the Q vector cal- 
culated on the proton and lead side, respectively. The lines 
show the result of combining the the two event planes. 
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FIG. 15. (Color online) Comparison of the elliptic flow co- 
efficients for |?7| < 1 calculated from the second cumulant 
method (solid line), from the scalar product method with 
event planes defined by charged particles in 2 < \r\\ < 2.5 
and 0.3 < p± < 3.0 GeV (dotted line), and from the scalar 
product method with event planes defined by the transverse 
energy in 3.2 < J77J < 4.8 (dashed line). 



r] limits defining the Q vectors. The flow coefficients are 
then calculated as 



.A.B 



{SP} 



< Qpf' cos {n{4>k 



*«))> 



Wn 



(16) 



The flovif coefficients with reduced statistical error can be 
obtained with combined event planes on the proton and 
Pb sides. 

In Fig. [13] we show the elliptic and triangular flow co- 
efficients obtained from the scalar-product method with 
the Q vector from the bin 2.0 < |?7| < 2.5 on either 
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FIG. 16. (Color online) Same as Fig. [15] for the triangular 
flow. 



the proton or the lead side. We notice that the two re- 
sults are very consistent, with slightly smaller statistical 
errors for the Q vector defined on the lead side. This re- 
flects a better resolution of the event plane in that case. 
The azimuthally asymmetric initial source for hydrody- 
namic evolution is longitudinally extended, which yields 
a strong correlation between the event planes on the lead 
and proton sides. The observed two-particle correlation 
functions are almost symmetric for A77 > and A77 < 0, 
which shows that the correlations are similar on the pro- 
ton and the lead side [4| . The consistency of flow correla- 
tions deflned with Q vectors for positive and negative ra- 
pidities justifies the use of the combined Q vector, which 
reduces the statistical error. 

The results obtained with the Q vectors defined by 
charged particle tracks or the calorimeter energy are com- 
pared to the results of the second cumulant method in 
Figs. [15] and [161 The elliptic and triangular flow coef- 
flcients obtained from the different deflnitions of the Q 
vector are very similar. The second cumulant harmonic 
flow is calculated for smaller average pseudorapidity sep- 
aration of the pair, thus contains some contribution of 
non-flow effects which increase the observed correlations. 
We expect that in the presence of additional non-flow 
correlation in the small system, such deviations could be 
larger. By comparing the second cumulant w„ to meth- 
ods using large rapidity gaps, the importance of such 
non-flow correlations could be estimated in the data. 



C. Correlations of flow and the initial geometry 

One of the main reasons to study the flow is to acquire 
the knowledge on the early phase of the reaction. One 
result that holds event-by-event is the proportionality of 
the eccentricity coefflcients of the "geometric" distribu- 
tion, £„, to the coefficient of the harmonic flow of the 
produced hadrons, «„. In Figs. [T7] and [TH] we show the 
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FIG. 17. The scattered plot of the event-by-event eccentricity- 
ehiptic flow correlations. Glauber+NB, correlation coefficient 
0.85. 
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FIG. 18. Same as Fig. [T7]for the triangular case. Correlation 
coefficient 0.74. 



event-by-event scattered plots of eccentricity-flow distri- 
butions. For this calculation the hydrodynamic events 
are combined from 1000 THERMINATOR events corre- 
sponding to the same freeze-out configuration. We notice 
large correlation coefficients, defined as 



P 



{tnVnj'i.}) - {en){Vn{'2.}) 

var(e„)var(w„{2}) 



(17) 



in these distributions, p — 0.85 for the elliptic and p ~ 
0.74 for the triangular case, respectively. This feature, 
well know for the A-A collisions, is vividly present in our 
treatment of the p-Pb collisions. 



V. CONCLUSION 

We have analyzed various aspects of soft collective 
dynamics of the relativistic p-Pb collisions in the ap- 
proach consisting of three stages: Glauber modeling of 
the initial phase, event- by- event viscous 3-1-1 D hydrody- 
namic, and statistical hadronization. Our analysis shows 
that the collective dynamics may very well be present 
in the highest-centrality p-Pb system formed in ultra- 
relativistic heavy-ion collisions. The application of the 
three-stage model, where the shape fluctuations in the 
initial stage are carried over to the harmonic flow coeffi- 
cients in the hadronic spectra, allows for a quantitative 
understanding of the data for V2 and V3, as well as to de- 
scribe the ridge structures in the two-particle correlation 
functions. The issues connected to the femtoscopic vari- 
ables in p-Pb collisions, which display considerable sensi- 
tivity to collectivity, have been presented elsewhere [24| . 

Thus, following the successful experience of describ- 
ing the A-A collisions in the three-stage approach, we 
argue that the collective scheme provides a natural and 
conventional explanation of numerous aspects of the soft 
dynamics of the "small" p-Pb system. 

In central p-Pb collisions, the initial shape eccentric- 
ity parameters e„ are entirely due to fluctuations. These 
fluctuations are enhanced by the distribution of overlaid 
weights on the spatial distribution of the participant nu- 
cleons. We find it quite remarkable that the same distri- 
bution that explains the multiplicity distribution of the 
produced hadrons in minimum-bias collisions leads also 
to quantitative agreement for the values of the elliptic 
and triangular flow coefficients measured recently by the 
ATLAS collaboration ^]. This agreement includes the 
elliptic flow coefficient obtained from the four-particle cu- 
mulants. 

We argue that the lowest harmonic flow coefficient 
may be measured directly in the LHC p-Pb experiments, 
hence we compute them through various methods (cumu- 
lant, scalar-product, rapidity-gap). These predictions, as 
well as the femtoscopic radii ^], will hopefully be veri- 
fied shortly in the upcoming experimental analyzes. 
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